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I . FORMULATION 

The transition to turbulence in boundary layers is investigated by direct numeri- 
cal solution of the nonlinear, three-dimensional, incompressible Navier-Stokes equa- 
tions in the half-infinite domain over a flat plate- Periodicity is imposed in the 
x-direction (strearawise) and in the z-direction (spanwise) ; the y-coordinate extends 
from 0 to ». 

The simulation is periodic in the x-direction, unlike the experiments in which 
the boundary-layer thickness grows in the x-direction* A forcing term is added to 
the x-momentum equation that approximates the convection terms associated with this 
spatial growth. The approximation is based on boundary- layer assumptions (applied 
only to the mean flow) and the self-similarity of the mean-velocity profile* With 
this forcing applied, the laminar velocity profile, instead of becoming an error 
function and thickening without bounds, is a Blasius profile- Thus, the stability 
characteristics are very close to the experimental characteristics. Furthermore, the 
equation can be written in a moving reference frame, so that the boundary- layer 
thickens in time while retaining a Blasius profile* The procedure of adding a forc- 
ing term allows the disturbances to extract energy from the mean flow, and is much 
preferable to a procedure in which the mean-velocity profile would be imposed. 

The spatial representation is spectral in all directions [1] . The basis func- 
tions that are used represent divergence-free velocity fields and satisfy the bound- 
ary conditions as suggested by Leonard and Wray [2]. Leonard and Wray applied a weak 
formulation, which eliminates the pressure and allows an accurate and straightforward 
time-advance scheme. Leray f s weak formulation is used here [3]. An advantage of 
this formulation over Leonard’s is that it keeps the numerical Stokes operator real, 
symmetric, and negative-definite. On the other hand, Chebyshev polynomials cannot be 
used . 

The x- and z-directions are treated by Fourier series. In the y-direction, 
the velocity field is first split into "irrotational" and M vortical M components in 
order to better accommodate two different length scales. The length scale of the 
vortical component is 6, the thickness of the boundary layer. The length scale of 
the irrotational components is A, the wavelength in the (x,z) plane, which is sig- 
nificantly larger than 6. This irrotational component can be represented by a 
single exponential function for each horizontal wave-vector. To represent the vorti- 
cal component, an exponential mapping is applied from [0,®[ into [0,1[, and shifted 
Jacobi polynomials are used in the transformed coordinate. The vortical component is 
infinitely differentiable over the closed interval [0,1], so that the convergence of 
the polynomial method will be faster than algebraic. The cost of the transforms from 



real space to Jacobi space is of the order of N 2 . Figure 1 is a plot of the first 
few basis-functions versus y. All the functions decay exponentially as y « but 
the first function, which includes the irrotational component, decays much more 
slowly than the other ones. 

The time-advance scheme is hybrid and second-order accurate. The convection 
terms are treated by a Runge-Kutta scheme which is explicit, third-order accurate, 
and conditionally stable; the Stokes terms are treated by the Crank-Nicolson scheme. 

II. RESULTS 

In order to check the convergence of the method, the Orr-Sommerfeld equation 
was solved for a Blasius profile and for a real wave-number. This problem is known 
to produce a few discrete eigenvalues and a continuous spectrum on the C r * 1, < 0 

axis [4]. Figure 2 is a contour plot of the error in the principal discrete eigen- 
value as a function of Yo, the length scale of mapping, and Ny, the number of points 
in the y-direction. The convergence as Ny -*■ » with Yo fixed is very fast. It is 
expected to be faster than algebraic, but not as fast as exponential [5]. The plot 
also indicates the optimum value of Yo: about 26*. Figure 3 shows that the numeri- 

cal spectrum includes a string of eigenvalues that becomes denser and tends to the 
C r * 1, < 0 axis. Its convergence is much slower than that for the discrete 

eigenvalues; the reason is that the corresponding eigenfunctions behave like sine 
waves as y 00 , which makes them hard to approximate with the expansion functions 
in Fig. 1. 

The early nonlinear stages of transition of a Blasius boundary layer, disturbed 
by a vibrating ribbon, were then simulated m three dimensions. A two-dimensional 
Tollmein-Schlichting (TS) wave of finite amplitude was introduced m the initial 
field, as well as three-dimensional white noise of much lower energy. The streamwise 
period was twice that of the TS wave; the spanwise period was chosen much longer to 
avoid constraining the spectrum. Spanwise lines of particles were introduced, near 
the critical layer, to simulate the smoke lines used m experiments. 

Figure 4 summarizes the time-evolution of the flow. The energy of the funda- 
mental TS wave and the energy carried by all the other wave-vectors are plotted sep- 
arately, The TS wave grows from branch I to branch II of the TS stability diagram, 
then starts decaying. The energy of the other modes remains small until after the 
flow crosses branch II; then It grows very rapidly, and nonlinear interactions take 
place. The shape factor H of the boundary layer remains at the Blasius value of 
2.6 until transition occurs; then it rapidly decreases. The agreement with Kachanov ? s 
experiments is excellent [6]. 

Simulations were conducted with the same background noise, but different values 
for the TS wave amplitude. Figure 5 contains top views of the particles in the bound- 
ary layer. If the maximum TS wave amplitude is less than 0.3%, transition does not 
occur. In Fig. 5(a), with an amplitude of 0.9%, three-dimensional breakdown occurs 
and is of the subharmonic or "H M type (the lambda-shaped particle lines are staggered) . 



In Fig. 5(b), with amplitude 5Z, the lambda patterns are not staggered, indicating a 
Klebanoff-type breakdown. The patterns appear "broken,” a result of the randomness 
of the initial three-dimensional disturbance. The qualitative agreement with Saric's 
experiments is good [7]. 

Figure 6 is a plot of the spectrum in an (x,z) plane at the beginning of an 
H-type breakdown. The fundamental TS wave still dominates the spectrum; its higher 
harmonic is also present. The growing three-dimensional subharmonic component is 
obvious; the wave number and the broadband character of the instability agree very 
well with Herbert's small disturbance theory [8]. 
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Fig. 1 First four basis-functions. Fig. 2 Error in Orr-Sommerfeld 

eigenvalue. 
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Fig. 6 Two-dimensional spectrum. 
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